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ABSTRACT 

'  An  algorithm  for  structured,  large-scale,  convex  quadratic  programming 
problems  is  described.  The  structure  of  the  constraint  matrix  is  block 
diagonal  with  a  small  number  of  coupling  constraints  and  variables.  The 
algorithm  utilizes  twice  a  decomposition  procedure  that  was  developed  earlier. 
The  first  time  the  decomposition  procedure  is  used  to  break  up  the  coupling 
constraints,  and  the  second  time  it  is  used  to  break  up  the  coupling  variables. 
Preliminary  computational  results  are  also  reported. 

(•  /  ) 
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SIGNIFICANCE  AND  EXPLANATION 


•The  block  diagonal  structure  with  a  small  number  of  coupling  constraints 
and  variables  usually  arises  from  the  formulation  of  multitime  period  and  multi¬ 
division  production  scheduling  and  distribution  models  in  large  corporations. 
Each  block  is  concerned  with  the  operation  of  one  division  during  one 
time  period  considered  in  isolation.  The  coupling  constraints  arise  from  the 
use  of  common  resources  of  all  divisions  or  from  combining  the  output  of 
divisions  to  meet  overall  demands.  The  coupling  variables  represent  acti¬ 
vities  that  affect  the  operation  of  divisions  in  more  than  one  time  period. > 

In  this  paper  we  propose  a  decomposition  method  for  convex  quadratic 
programming  problems  having  block  diagonal  structure  with  a  small  number 
of  coupling  constraints  and  variables.  The  basic  idea  of  decomposition  methods 
is  to  break  up  a  large  complex  problem  into  a  sequence  of  subproblems  that 
are  easier  to  solve* 


AN  ALGORITHM  FOR  STRUCTURED,  LARGE-SCALE  QUADRATIC 
PROGRAMMING  PROBLEMS 


Cu  Duong  Ha 


1.  Introduction 


The  problem  that  we  shall  be  concerned  with  is 


r  . 

1  (*e.,x.>  ♦-<*.  ,C.x.  >) 

i-0  11  Z  i  i  i 


subject  to  BjXq  +  DjXj 


♦  D.x. 
2  2 


Vo  +  y  p  -  dp 

Vo  +  A1X1  +  A2x2+-  *  *+  ApXP  "  * 

x.  >  9  i  -  0,1,. ..,p 

“i  «•  m.  m 

where  x£  £  K  ,  c.  €  K  ,  di  e  *  1  ,  «  c  *  0  ,  A.,  C..  and 

D.  are  matrices  of  dimensions  l  x  n.,  m.  x  nQ,  n.  x  n.  and  m.  x  n. 

respectively;  and  <•,*>  is  the  inner  product. 

P 

We  denote  x:»  (x^,Xj,...,x  )  ,  n;»  J  n. 

P  i=0  1 

We  assume  that  the  problem  is  feasible  and  is  symmetric  and 
positive  semidefinite  for  all  i. 

Note  that  we  do  not  rule  out  the  case  C^=0  for  all  i  ;  accordingly 

our  algorithm  can  be  used  to  solve  block  diagonal  linear  programming 

problems  with  coupling  constraints  and  variables. 

It  is  well-known  that  decomposition  approaches  for  large-scale 

problems  which  uses  dual  methods  have  several  disadvantages  (see,  for 

example,  [Shapiro,  1979]  or  [Lasdon, 1970] ) .  To  circumvent  these 

drawbacks  we  developed  a  decomposition  procedure  by  combining  a  dual 
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method  and  the  proximal  point  algorithm  and  used  it  for  solving  block 
angular  linear  programming  problems  with  coupling  constraints  and  a 
large-scale,  nonlinear  problem  of  structural  engineering  ([Ha,  1981]  and 
[Kaneko  and  Ha,  1980]).  In  this  paper  we  apply  the  same  decomposition 
procedure  to  the  problem  (1.1).  Ihe  resulting  problem  is  a  quadratic 
programming  problem  having  only  coupling  variables;  and  its  dual  problem 
has  only  coupling  constraints.  Then  we  modify  our  decomposition  procedure 
to  break  up  that  dual  problem  completely  into  small  subproblems. 

The  organization  of  this  paper  is  as  follows.  Section  2  is  the 
review  of  the  decomposition  procedure  in  [Ha,  1981] .  Section  3  is  the 
proposed  solution  method  for  the  problem  (1.1).  Section  4  is  the  imple¬ 
mentation  and  computational  results  of  the  algorithm. 


decompostion  procedure 


1981] ) 


In  this  section  we  consider  tile  convex  programming  problem 


min  f(x) 

subject  to  (2.1) 

X  £  c 
Ax  =  a 

n  m 

where  x  gIR  ,  a  €  E*  ,  A  is  an  m  x  n  matrix,  f  is  a  convex  function  defined 
on  IRn,  and  C  is  a  nonempty  closed  convex  set  in  IRn.  We  assume  that  the 


problem  (2.1)  has  optimal  solutions 


3 


Let  {A^)  be  a  sequence  of  positive  numbers  bounded  away  from 

zero  (i.e.jX^  >_  A  >  0  for  all  k).  Let  x°  be  an  arbitrary  point 

n  _  _  ,  r  0  1  k-,  k  +  1  , 

in  TR  .  Suppose  we  already  have  ix  ,  x  ,  ...  ,  x  j  ,  then  x  is 

the  unique  optimal  solution  of  the  problem 

min  f(x)+  -ry-  [j  x  -  x*  ||  2 
k 

subject  to  (2.2) 

x  e  c 

Ax  =  a 

Then  by  the  proximal  point  method  ( [Rockafellar,  1976  a  and  b]  and 
[Bre*zis  and  Lions,  1978])  the  sequence  {x  }  will  converge  to  an  optimal 
solution  of  the  problem  (2.1) . 

The  advantage  of  (2.2)  over  (2.1)  is  that  the  objective  function  of 
(2.2)  is  strongly  convex  although  the  function  f  may  not  be  so.  A 
function  h  defined  on  a  convex  set  S  is  said  to  be  strongly  convex 
with  modulus  a  >  0  if 

h((l-t)x  +  tz)  <  (l-t)h(x)  +  t  h(z)  -  i  a(l-t)t  ||  x  -  z  ||2 

for  all  x  and  z  in  S  and  0  <  t  <  1 . 

In  order  to  solve  the  problem  (2.2)  we  use  a  dual  method.  The 
dual  of  (2.2)  is 


max  g(y) 

y  € 


(2.3) 


! 


where 


g(y):=  min  (f(x)  +  [jx  -  xk|j2  -  <  y  ,  Ax  -  a>)  (2.4) 

x(?C  Alc 

1  k  2 

Since  f(x)  +  — —  ||  x  -  x  ||  is  strongly  convex.  The  problem  (2.4) 
k  m 

has  a  unique  optimal  solution  for  any  y  in  ®  .  In  other  words,  g(y) 
is  finite  everywhere.  It  can  also  be  proved  that  the  function  g  is 
Lipschitz  continuously  differentiable  and  the  derivative  g'  of  g  at  a 
point  y  is  given  by 

g' (y)  =  Ax(y)  -  a 

where  x(y)  is  the  optimal  solution  of  (2.4)  corresponding  to  the  given  y. 

If  y*  is  the  optimal  solution  for  (2.3)  then  x(y*)  is  the  optimal 
solution  for  (2.2). 

In  summary,  the  general  scheme  for  solving  (2.1)  is  to  choose  points 

x°  €  ®n,  y  €  and  then  solve  the  minimization  problem  (2.4)  to  obtain 

x(y) ,  g(y)  and  g' (y) .  If  y  is  not  an  optimal  solution  for  (2.3),  update  y 

and  solve  (2.4)  again  until  the  optimal  solution  y*  of  (2.3)  is  obtained. 

If  x(y*)  is  the  same  as  x®  then  it  is  an  optimal  solution  for  the  original 

problem  (2.1);  otherwise  replace  x°  by  x(y*)  and  repeat  the  procedure. 

Note  that  in  (2.4)  we  do  not  have  the  constraints  Ax  =  a  which  can  be 

the  coupling  constraints  in  the  case  that  (2.1)  is  a  structured  large- 

scale  problem.  Because  of  computational  experience  with  {A^}  in 

[Ha,  1981)  we  shall  keep  A  fixed  ; i.e. , A  =  X  >  0  for  all  k. 

k  k 


We  use  the  multilevel  concept  to  explain  our  method  of  solution 

for  the  problem  (1.1).  It  is  a  three  level  algorithm.  The  first  level 

is  the  application  of  the  decomposition  procedure  above  to  the  problem 

0  0  0  0. 

(1.1).  Starting  from  a  given  point  x  =  (x.  ,  x,  ,  ,  x  )  We  find 

0  1  p 

the  unique  minimizer  x  =  (xQ  ,  x^  ,  ...  x^)  for  the  following  problem 


I  (<V*i>  ♦  i  <VCi*i>  *  I 
1=0 


1  I  °ll2* 

2A  ilxrxill  > 


subject  to  BjXq  + 


B2X0  +  °2X2 


+  D  X 
P  P 


Vo  +  Aixi  +  * 


.  +  A  x 

P  P 


0  i  =  0,1,  . . . ,  p  . 

If  x  is  close  enough  to  x"  ,  x  is  approximately  an  optimal  solution 
for  the  original  problem  (1.1)  ?  if  not,  we  replace  x^  by  x  and  solve 

(3.1)  again.  We  solve  (3.1)  by  a  dual  method.  The  dual  problem  of 


(3.1)  is 


max  g(y) 


(3.2) 
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where 

g(y): 


Let 

and 

where 

g(y) 

here. 


min 


i«0 

subject  to  BjXq  +  D^x^ 
B2X0 


B  xft 
P  o 


f  «ci,x.> +i<xi,C.x.>  +  ^  |jx.-x°||2)+<y,(  l  A.x.-a) > 
i-0  1  x  *  1  1  1  ZA  1  1  i=0  1  1 


+  V2 


(3.3) 


+  D  x  ■  d 
P  P  P 


*•  >;  0  for  i  =*  0,1,..., p  . 

-  1  0  T 

c.:“  c.  -  yx-  +  A.y 

l  l  A  i  i7 

C.:»  C.  +  i  I 

l  l  A 

I  is  the  identity  matrix  with  appropriate  dimension;  then 
can  be  rewritten  as 

g(y)--<..y>*i||x?||2.8l(y>  , 
gj(y)  is  the  optimal  objective  value  of  the  problem: 


min 

l  (<c.,x. > 
i=0  1  1 

+  2  Xi'Vi> 

subject  to 

Vo  +  Vl 

'  dl 

B2x0 

• 

+  °2X2  *  d2 
•  • 

• 

B  xA 

P  o 

•  • 

♦  D  x  *  d 

P  P  P 

x.  >  0 

i  ” 

for  i  *  0,1, ... ,p  . 

(3.4) 


For  a  given  y,  the  second  level  will  give  us  x(y)  the  optimal 
solution  of  (3.4).  Then  we  can  compute  g(y)  and  g'(y).  If  y  is  the 


optimal  solution  of  (3.2)  then  x(y)  is  the  variable  x  that  we  are  looking 
for;  otherwise,  we  can  use  any  unconstrained  optimization  algorithm  to 
update  y  and  send  the  new  value  of  y  to  the  second  level.  The  second  level 
is  concerned  with  the  problem  (3.4)  for  a  given  y  receiving  from  the 
first  level  (  c  is  a  function  of  y  ) .  The  ordinary  dual  of  (3.4)  •  . 

can  be  written  as  { [Mangasarian,  1969]). 


r  1  -  K 

min  I  Or<x.,C.x.»  +  £<v.,d.> 

(x,v)  i=0  1  111  i=l  1  1 


(3.5.1) 


subject  to  +  D*v^ 


C2X2  +  D2v2 


C0X0  +  Blvl  + 


>  -c. 


*-C2 


C  x  +  D  v  >  -c 
P  P  PP“  P 
T 

+  B  v  >  -cn 

p  p  -  0 


(3.5.2) 


(3.5.3) 


The  problem  above  is  a  quadratic  programming  problem  having  block 
angular  structure  with  (3.5.3)  considered  as  the  coupling  constraints. 

It  can  be  solved  by  the  decomposition  procedure  outlined  in  Section  2. 
But  in  this  special  case,  because  the  objective  function  (3.5.1)  is 
already  strongly  convex  with  respect  to  x,  we  can  simplify  the  procedure 
by  adding  only  quadratic  terms  for  the  variable  v.  More  specifically, 
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let  {vk>  be  a  sequence  of  positive  numbers  bounded  away  from  zero  and 

0  0 
v  be  a  starting  point.  Suppose  we  have  generated  from  v  the  sequence 

.  1  1  .  2  2  k  k.i  ,  k+1  k+1  . .  .  ,  . 

{  (x  ,v  ) , <x  ,v  ) ,  . ..,  (x  ,v  )j  ;  let  (x  ,  v  )  be  the  optimal  solution 

of  the  following  problem 

P 


min 

(x,v) 


I  (|<x.,Cx  >)+  l  <vi,di>  +  2J7*  ||  v-vk  || 2 
i=0  Z  1  1  1  i-1  1  1 


subject  to  CjXj  + 


C2X2  +  D2V2 


>  -c,  (3.6) 


Vo  +  Bivi  + 


C  x  +  D  v  >  -c 
P  P  P  P"  P 

+  BTv  >  -c 
P  P  "  o 


k  k 

then  the  sequence  { (x  ,  v  )  }  will  converge  to  an  optimal  solution  for 
(3.5)  (justifications  will  be  given  at  the  end  of  this  section).  Again, 
we  solve  (3.6)  by  a  dual  method.  The  dual  of  (3.6)  is 


max  h(u) 
u>0 


(3.7) 


where,  for  a  fixed  u  , 

P 


h(u)  =  min 
(x,v) 


+  <u’"cO~Vo“  *  Vi> 

j=l  J 


subject  to  C^x^  + 


C2X2  +  °2V2 


>-i 

> 0.8) 

-  2 


C  x  +  D  v  >  -c  . 
P  P  P  P”  P 
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The  third  level  is  to  solve  the  minimization  problem  (3.8)  above  for 
given  u  and  v  and  then  pass  back  to  the  second  level  its  optimal 
solution  (x(u)  ,  v(u)).  The  problem  (3.8)  can  be  separated  into  (p+1) 
small  subproblems.  The  subproblem  corresponding  to  i=0  is 


minl(xo*Vo>  "  <u‘Vo>  * 


(3.9; 


whose  solution  is  =  u;  for  i  >  1  ,  it  is 

(*miv  )  I<VVi>+<Vi>  +  2^llvi“Vi112  ■  <u>BiV  (3. 


10) 


l  l 


subject  to  C.x.  +  D.v.  >  -  c. 

J  li  ii“  l 

which  can  be  solved  by  any  algorithm  for  quadratic  programming  problems. 

In  the  remainder  of  this  section  we  shall  prove  a  result  to  justify 
the  second  level.  We  consider  the  problem 


min  9(x)  +  <J>(v)  (3.11) 

(x,vKC 

where  C  is  a  nonempty  closed  convex  set,  0  and  <p  are  closed  covex 
functions  having  real  values. 

We  suppose  the  problem  has  optimal  solutions. 

Lemma  3.1 

0  1 

Let  y  be  a  fixed  positive  number.  Given  any  two  points  v  and  v  , 


U°.V0):. 


argroin  9(x)  ♦  4>(v)  ♦  i  ||v-V°||2 
(x,v)cC  Zy 


we  define 


then 


(x  ,V  ):*  argmin  0(x)  +  4>(v)  +  -577  !) v— V  ||  , 


(x,v)cC 


2y 


||v°-vl||2  >  ||  v0-^1 1| 2  +  ||v°-v°+v1-v1||2 


(3.i: 


and 


11  v°— v1 1|  >  llv0-^1  II 


(3.i: 


Proof 

It  is  clear  that  (3.13)  follows  from  (3.12),  so  only  (3.12) 
needs  to  be  proven.  We  have 

II  ^-v1 1| 2  =  ||v0-v°+v°-v1+vl-v1|}2 

*  II  v0-^1 1| 2  +  ||v°-v°+v1-v1||2  +  2(v°-v1,v0-v°+vl-v1>  , 

so  to  prove  (3.12)  we  shall  show  that 

<v°-vl,  v°-v°  +  v1-v1>  >  0  . 


Let 


F(x,v):=  0(x)  +  $(v)  +  ^(x,v) 
and  G(x,v):“  (2y)  *  |Jv-V^||2  . 


Since  F  and  G  are  both  proper  convex  and  the  effective  domain  of 
G  is  the  entire  (x,v)-space,  we  have  ([Rockafellar,  1968]) 


It  follows  from  the  definition  of  the  subdifferential  that 


(0,0)  e  3(F  +  G)  (i°,V°) 

*  3F(x°,V°)  +  3G(x°,V°) 

-  3F(i°,V°)  +  {0}  x  p*1(v°-V°) 
or 

(0>y”1(V°-V0))  e  3F(x°,V°)  . 

Similarly,  we  have 

(0,p"’1(V1-V^))  €  BFG^.V1)  . 

Therefore  by  the  monotonicity  of  3F 

(o-o,x°-xl>  +  y'1  < v^v^v^v1,^0-^1)  >0  , 
or  <  V^V1 ,  V^V^V^V1  >  >  0  . 

Theorem  3.1 

In  addition  to  the  above  assumptions,  we  assume  that  6  is 
strongly  convex.  Let  (y^)  be  a  sequence  of  positive  numbers  bounded 
away  from  zero  and  let  v^  be  an  arbitrary  point.  The  sequence 
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(xk+*,vk+*):»  argmin  0(x)  +  4>(v)  +  -=7— ||v-vk||  2 
(x,v)eC  2pk 

converges  to  a  point  (x*,v*)  which  is  an  optimal  solution  for 
<3.11). 


Proof 

Let  (x,v)  be  an  optimal  solution  of  (3.11) .  It  is  easy  to  see 
that  (x,v)  is  also  the  optimal  solution  of  the  following  problem 


min 

(x,v)eC 


0(x)  +  4>(v)  +  2^||v-v||2 


Applying  lemma  3.1  with  y  =  y^  .  =  v  ,  s  v*1  and 

-1  k+1 

V  *  v  ,  we  have 

||vk+*-v|j  <  ||vk-vj|  (3.14) 

and  l|vk+1-v||2  +  ||vk+1-vk||2  <  ||vk-v||2  . 


k 

The  first  inequality  gives  us  the  boundedness  of  the  sequence  {v  }  . 
The  second  shows  that 


I  l|vk,l-vk||2 


k=l 


<  00  . 


It  implies  that  ||vk+^-vkJ|  ■*  0 
to  show  that  the  sequence  {x  } 


as  k  -*•  00  .  Now,  we  are  going 
is  also  bounded.  First  we  have 


a/  k+l\  At  k+1*  ^  at  k+l»  ..  k+lv  ,  1  II  k+1  k*  .,  k. 

8(x  )  ♦  ♦(v  )^0(x  )  +  ♦(v  )+r— -||v  II  ^  6(*  )  +  $(v  )  , 

ZMk 

so  8(xk+*)^  0(x*)  ♦  $(v*)  -  $(vk+l)  . 

k.  1c 

The  boundedness  of  {v  }  implies  the  boundedness  of  {$(v  )}  •  bet 

M  be  a  constant  such  that 

♦(v1)  -  $(vk)  <  H  Vfc  , 

then 


xk«W:*  {x|  8(x)  <  8(x*)  +  M}  . 


The  set  W  is  compact  by  the  strong  convexity  of  l|>  ;  consequently 
{x^}  is  bounded. 

k  k 

Since  the  sequence  {(x  ,v  )}  is  bounded,  there  is  a  subsequence 
k.  k. 

{(x  J,v  J)}  that  converges  to  a  point  (x*,v*)  e  C  .  For  any  point 
(x,v)  e  C  ,  by  lemma  3.1  we  have 


k  *  k . 

0(x)  +  4>(v)  >  0(x  +  <J>( V 


.  k .  -1  k . 

+  — <v  J  -v  J,  v 

Tc. 

J 


Letting 


k. 

J 


we  have 

k.  k. 

0(x  ^)  ♦  $(v  -*■  0(x*)  +  4>(v*) 

k.-l  k. 
y  J  -  y  ^  +  o  . 


Since  {u.  }  is  bounded  away  from  zero 


f 


I 


.  k.-l  k.  k. 

— — (  v  ^  -v  ^,v-v  ^  )  -*•  0  , 

Tt. 

J 


so 


e(x)  +  <j> c v)  >  0(x*)  +  4»Cv*)  . 


This  is  true  for  an  arbitrary  (x,v)  t  C  ,  so  (x*,v*)  solves  the 
problem  (3.11).. 

Using  (3.14)  with  v  replaced  by  v*  we  have 

||vk+1-v*||  <  ||vk-v*||  , 

so  the  sequence  {v  }  has  to  converge  to  v*  . 

Suppose  {x  }  has  two  cluster  points  x  and  x  ,  then  both 
(x,v)  and  (x,v)  are  optimal  solutions  for  the  problem 

min  0(x)  +  $(v)  +  ^-||v-v||2 
(x,v)cC  Zy 


1  -  2 

for  any  positive  number  y  .  But  the  function  8(x)  +  $(v)  +  || v-v|| 

is  strongly  convex  with  respect  to  (x,v)  ,  so  the  problem  above  has 
a  unique  optimal  solution.  Hence  x  =  x  and  the  sequence  {(x  ,v  )} 
converges  to  (x*,v*)  . 


4.  Implementation  and  computational  results 

At  Level  1  the  problem  to  be  solved  is  (3.2)  which  is  a  nonlinear, 
unconstrained  maximization  problem,  whose  objective  function  g(y)  is 
differentiable  and  its  derivative  g* (y)  is  available.  To  solve  (3.2) 


\ 


IS 


we  use  the  Broyden-Fletcher-Golfard-Shanno  (BFGS)  algorithm,  which  is 
available  in  the  Harwell  Subroutine  Library  and  known  as  VA13A  ( [MACC,  1976] ) . 
For  Level  2  we  have  to  handle  the  nonlinear  maximization  problem  (3.7) 
with  nonnegativity  constraints.  In  order  to  use  an  unconstrained 
algorithm  we  modify  (3.5),  (3.6)  and  (3.7)  as  follows. 

We  convert  the  inequality  coupling  constraints  (3.5.3)  into 
equality  constraints  by  adding  slack  varibles  w;  (3.5)  then  becomes 


P  .  P 

min  l  (i<x.,C.x.  >)  +  J<v.,d.  > 

(x,v),v  i=0  2  1  1  i=l 


subject  to  +  Djv2 


>  -c. 


c2*2  +  V2 


-  ~c2  (3.5*) 


C  x  +  D  v  >  -c 

P  P  P  P  ”  P 


Vo  +  B1V1  + 


♦  B  v  -  w  »  -cn 
P  P  0 


v  >  0 


Instead  of  (3.6)  and  (3.7)  we  have 


r 
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Its  ordinary  dual  problem  is 


mflir 

(Xj.vJ.Z 


sT_ 


subject  to  C^x^  -  C^z  =  0 


(4.2) 


-f-  v.  -  D.z  =  -d. 

\  1  1 

z  >  0 


From  the  constraints  we  have 


-  -T 

C.x.  “  C.z 
ii  i 


Since 


is  symmetric  and  positive  definite,  this  implies 


x.  ■  t  . 
i 


We  also  have 


v^  =  ^(D^z  -  d^) 

*  Wi '  3i) 

Replacing  z  by  x^  and  v^  by  the  expression  above  in  (4.2)  , 
(4.2)  is  equivalent  to 


max 

x.>0 

i" 


-  4< 


|<xi,Cixi>  -  ~<Vk(Dixi“d.),  ^(DjX.-dJ)  -(cj.x,^) 
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or 


min 

x.>0 

i"* 


(C.+y,  dTd.)x.)  +  (x., 

i  k  1  i  i  i 


(4.3) 


Remark 


1  }c  2 

If  at  Level  2  we  also  add  the  quadratic  term  [|x-x  ||  into 

ZMk 

the  objective  function  (using  the  proximal  point  algorithm);  the  sub¬ 
problem  (3.10)  of  level  3  would  be: 

min  l<xi*^ixi>  +<vi'di>  +  2uT  llVi"ViH2+  2SrIlxi"xiH2-<u»Bivi> 

(x.,v.)  *  Tt 


i  i 


subject  to  +  ®ivi  £  ~c£ 


Following  the  same  arguments  as  above  one  still  could  convert  this 
problem  into  a  problem  whose  only  variables  are  x^  .  But  instead  of 
(4.3)  one  would  have 


min 

x. 

l 


T<*i-<V\ll*(I%lai1>DiDi<I*\lai1))*t> 

-lr-l^T,  -U  rlk  - 


subject  to  (b.  I+C.l)x.  >  C. 4x.k 
J  k  x  i  *•  l  i 


srl. 


which  is  much  more  complicated. 

There  are  several  algorithms  for  solving  the  quadratic  programming 
problems  (4.3).  With  pivoting-type  algorithms  we  cannot  use  the  optimal 
solution  of  the  previous  iteration  as  the  starting  point  for  the  current 
iteration.  The  Best-Ritter  algorithm  ([Best-Ritter,  19761)  which  is  a 


conjugate  direction  method  allows  us  to  do  that,  so  we  use  it  to  solve 
(4.3). 


2C 


As  in  the  case  of  A  we  keep  fixed?  i.e.,  \  =  W>0  for  all  k. 

In  that  case  the  matrices 


C.  +  V  D.  D. 
1  11 


(4.4) 


remain  unchanged  throughout  the  process,  so  they  need  to  be  computed  only 
once. 

We  wrote  a  computer  program  to  test  the  algorithm,  using  FORTRAN  V 
on  the  UNIVAC  1110  of  the  Madison  Academic  Confuting  Center  of  the 
University  of  Wisconsin-Madison.  In  the  program  we  set  the  stopping 
criteria  as  follows. 


then  x  is  considered  to  be  an  optimal  solution  for  the  problem  (1.1). 

e2  =  10  This  is  the  tolerance  of  the  BFGS  algorithm  (problem  2.2). 

That  means  a  solution  y  is  accepted  as  optimal  if  a  relative  change  of 
size  e2  in  the  components  of  y  does  not  reduce  the  objective  value. 

We  have  the  same  tolerances  for  Level  2. 

For  Level  3  we  set  ■  10 

We  generated  test  problems  by  predetermining  the  size  of  the 
problem,  the  size  of  blocks,  and  the  number  of  coupling  constraints  and 
variables;  then  generating  randomly  data  of  the  problem.  The  matrices  Aj, 


» .  r .» •  !■> 
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and  are  90%  dense.  Each  entry  of  those  matrices  is  a  pseudo¬ 
random  number  in  the  range  [-50,  50],  obtained  by  the  random  number 

routines  of  the  Madison  Academic  Computing  Center  [MACC,  1978].  Each 

T 

positive  semi -definite  matrix  c.  is  of  the  form  E.E.  ,  where  E.  is  a 

1  xxi 

randomly  generated  matrix.  The  cost  coefficients  are  pseudo-random 
numbers  in  the  range  [-10,  10].  To  be  sure  that  the  problem  is  feasible 
we  randomly  generated  a  sequence  of  integers  in  between  0  and  5  (considered 
as  a  feasible  point)  and  then  multiplied  them  in  a  proper  way  with  the 
matrices  A^,B^  and  to  obtain  the  right-hand  side  coefficients.  The 
structures  of  the  test  problems  are  given  in  Table  4.1  below. 


'problen\. 

Problem 

Size 

Number  Of 

Blocks 

Coupling 

Variables 

Coupling 

Constraints 

I 

50  x  100 

10 

5 

5 

II 

59  x  125 

16 

2 

3 

III 

70  x  150 

20 

2 

3 

Table  4.1  Structures  of  Test  Problems 

We  compared  the  computational  results  with  the  results  obtained  by 
using  the  LCPL  package  available  at  the  Madison  Academic  Computing 
Center  of  the  University  of  Wisconsin-Madison.  LCPL  is  a  program  for 
solving  linear  complementarity  problems  and  quadratic  programming  prob¬ 
lems.  It  was  developed  by  Tomlin  [1976]  at  the  Systems  Optimization 
Laboratory  of  Stanford  University.  It  uses  a  variant  of  Lemke's  method 
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such  that  only  the  nonzero  coefficents  are  stored  in  packed  form  and 
the  basis  inverse  is  maintained  in  standard  product  form.  Computa¬ 
tional  results  are  shown  in  Table  4.2. 


^^Proble^N^^ 

LCPL 

Our  Algorithm 

A=  10,  y  =  50 

A  =10,  y  =  100 

A=  10,  y  =  200 

1 

10. 

40. 

18. 

19. 

II 

14. 

33. 

36. 

(1) 

III 

(2) 

158. 

34. 

29. 

Table  4.2  CPU  Time  (secs)  of  Test  Problems 

(1)  Numerical  errors. 

(2)  The  run  stopped  after  255  iterations  because  of  the  insufficiency 
of  space  to  store  the  inverse  of  the  basis  in  product  form  (ex¬ 
ceeded  the  limit  of  8000  words  of  the  eta  file) . 

From  Table  4.2,  we  see  that  the  computation  time  of  the  de¬ 
composition  procedure  is  not  yet  as  good  as  that  of  LCPL.  Several  points 
are  relevant  to  this  comparison.  First,  our  computer  program  was 
written  for  the  purpose  of  solving  a  limited  number  of  test  problems 
and  investigating  some  of  the  conputational  aspects  of  the  method? 

consequently,  it  was  not  written  with  a  great  a  programming  efficiency  as  was 
LCPL.  Secondly,  the  decomposition  approaches  are  intended  to  handle  large 

problems?  however  our  test  problems  are  by  no  means  large. 


In  terms  of  memory  storage  requirements,  our  method  is  much  better 
than  LCPL.  The  requirements  of  our  method  are  as  follows  (we  do  not 
count  the  storage  for  problem  data  since  it  is  approximately  the  same 
for  both  methods) : 


Storage  requirements  (neglecting  small  items) 

Level  1:  We  need  to  store  x^  ,  x  and  five  other  vectors  of  the  same 

dimensions  n,  so  the  total  number  of  words  is  7n.  The  require¬ 
ments  for  a  BFGS  routine  with  the  dimension  of  the  variables 
being  i  are  3 i  +  i (X+  13)/2. 

^  _  jc+2  fc+i 

Level  2;  Storage  for  the  vectors  (v  ,w  ) , (v,w)  and  (v  ,w  )  is 


r 

l((£  m.)-*n  );  working  space  needed  is  2((£  mj+n0).  There  is  another 
i=0  1  i=l  1 


BFGS  routine  with  the  dimension  of  the  variables  being  nQ. 

Its  storage  requirements  are  3nQ  +  nQ(n0  +  13)/2. 

Level  3:  For  a  quadratic  programming  problem  with  the  dimension  of  the 

variables  being  N  and  with  only  nonnegativity  constraints,  the 

2 

storage  requirements  of  the  Best-Ritter  algorithm  are  2N  +  ION. 

To  accomodate  problems  of  sizes  n^,  i  =  1,  2,  p,  N  should  be 


(4.4).  They  are  symmetric,  so  the  requirements  are 


f  n. (n.+l) 

L  _ k _ 

i=0  2 


Storage  requirements  for  our  test  problems  are  shown  in  Table  4.3: 


^^Level'S^ 

Problem 

I 

II 

III 

1 

700 

875 

1050 

2 

60 

33 

33 

3 

250 

290 

345 

4 

60 

21 

21 

5 

2055 

2390 

2684 

Total 

2495 

3609 

4133 

Table  4.3  Storage  Requirements  (words)  for  Test  Problems 

From  Table  4.3,  we  see  that,  for  problem  III,  our  method  required  only 
4133  words  of  memory,  while  LCPL  used  up  the  limit  of  8000  words  for 
eta  file  alone  and  still  has  not  reached  the  optimal  solution  of  the 
problem  (see  Table  4.2),  With  larger  problems,  the  decomposition 
procedure  would  lead  to  even  greater  savings  in  memory. 
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